Spectral products are created using mathematical combinations of specific bands (wavelength components). These spectral products can be useful for identifying spatial variations in vegetation, water or urbanization. This notebook covers the following spectral products: Fractional Cover, NDBI, NDVI, NDWI, SAVI, EVI
import datacube
dc = datacube.Datacube(app='Spectral_Products')
import sys, os
os.environ['USE_PYGEOS'] = '0'
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import xarray as xr
from datacube.utils import masking
from dea_tools.plotting import rgb, display_map
from odc.algo import to_f32
### EASI tools
sys.path.append(os.path.expanduser('../scripts'))
from ceos_utils.data_cube_utilities.clean_mask import landsat_clean_mask_invalid, landsat_qa_clean_mask
from easi_tools import EasiDefaults
from easi_tools import notebook_utils
easi = EasiDefaults() # Get the default parameters for this system
Successfully found configuration for deployment "eail"
cluster, client = notebook_utils.initialize_dask(use_gateway=False)
display(cluster if cluster else client)
print(notebook_utils.localcluster_dashboard(client, server=easi.hub))
a2da3b55
| Dashboard: http://127.0.0.1:8787/status | Workers: 4 |
| Total threads: 8 | Total memory: 29.00 GiB |
| Status: running | Using processes: True |
Scheduler-b27d8ebe-785b-4f5d-b2c3-e99981e63dc1
| Comm: tcp://127.0.0.1:36521 | Workers: 4 |
| Dashboard: http://127.0.0.1:8787/status | Total threads: 8 |
| Started: Just now | Total memory: 29.00 GiB |
| Comm: tcp://127.0.0.1:36751 | Total threads: 2 |
| Dashboard: http://127.0.0.1:39309/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:40735 | |
| Local directory: /tmp/dask-worker-space/worker-227hjslp | |
| Comm: tcp://127.0.0.1:35963 | Total threads: 2 |
| Dashboard: http://127.0.0.1:42569/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:40263 | |
| Local directory: /tmp/dask-worker-space/worker-9raofrjh | |
| Comm: tcp://127.0.0.1:36497 | Total threads: 2 |
| Dashboard: http://127.0.0.1:42951/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:46399 | |
| Local directory: /tmp/dask-worker-space/worker-pfff9htp | |
| Comm: tcp://127.0.0.1:35369 | Total threads: 2 |
| Dashboard: http://127.0.0.1:34953/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:36405 | |
| Local directory: /tmp/dask-worker-space/worker-hrwxic0a | |
https://hub.eail.easi-eo.solutions/user/jhodge/proxy/8787/status
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
<botocore.credentials.Credentials at 0x7fd232798340>
# Select a Product
product = "landsat8_c2l2_sr"
# Select an analysis region (Lat-Lon) within the extents listed above.
# Select a time period (Min-Max) within the extents listed above (Year-Month-Day)
# This region and time period will be used for the cloud assessment
# Kigoma, Tanzania
# latitude = (-4.95, -4.84)
# longitude = (29.59, 29.75)
# Dar es Salaam, Tanzania
# latitude = (-7.0, -6.75)
# longitude = (39.1, 39.35)
# Mtera Reservoir, Tanzania
# latitude = (-7.22, -6.80)
# longitude = (35.60, 36.00)
# Mining Region near Obuasi, Ghana
# latitude = (6.0985, 6.2675)
# longitude = (-2.050, -1.8629)
# Efate Island, Vanuatu (near Port Vila)
# latitude = (-17.810, -17.691)
# longitude = (168.247, 168.349)
# Port Vila, Vanuatu
# latitude = (-17.8, -17.65)
# longitude = (168.25, 168.375)
# Solomons Island - Coconut Crop Damaged by Beetles - FULL
# latitude = (-9.4412, -9.4357)
# longitude = (160.0241, 160.0288)
# Solomons Island - Coconut Crop Damaged by Beetles - MIDDLE
# latitude = (-9.4391, -9.4372)
# longitude = (160.0257, 160.0278)
latitude = easi.latitude
longitude = easi.longitude
# Time Period (Year-Mon-Day)
time_extents = ('2015-1-1', '2015-12-31')
# The code below renders a map that can be used to view the region.
display_map(longitude,latitude)
/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
landsat_dataset = dc.load(latitude = latitude,
longitude = longitude,
time = time_extents,
product = product,
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
output_crs = 'EPSG:6933',
resolution = (-30,30),
dask_chunks = {'time':1}
)
# Displays the first few values of each data array to check the content
# Latitude and Longitude numbers = number of pixels in each dimension
# Time = number of time slices in the dataset
landsat_dataset
<xarray.Dataset>
Dimensions: (time: 36, y: 342, x: 322)
Coordinates:
* time (time) datetime64[ns] 2015-01-05T15:46:53.790263 ... 2015-12...
* y (y) float64 4.418e+06 4.418e+06 ... 4.408e+06 4.408e+06
* x (x) float64 -7.386e+06 -7.386e+06 ... -7.376e+06 -7.376e+06
spatial_ref int32 6933
Data variables:
red (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
green (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
blue (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
nir (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
swir1 (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
swir2 (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
pixel_qa (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
Attributes:
crs: EPSG:6933
grid_mapping: spatial_refcloud_mask = masking.make_mask(landsat_dataset['pixel_qa'], cloud='not_high_confidence', cloud_shadow='not_high_confidence')
land_mask = masking.make_mask(landsat_dataset['pixel_qa'], water='land_or_cloud')
# Land and Water Dataset = Land and Water pixels with NO Clouds and NO Cloud Shadows
land_and_water_dataset = landsat_dataset.where(cloud_mask)
# Land Dataset = Land ONLY pixels with NO Clouds, NO Cloud Shadows and NO Water pixels
land_dataset = landsat_dataset.where(land_mask)
from ceos_utils.data_cube_utilities.dc_mosaic import create_median_mosaic, create_max_ndvi_mosaic, create_hdmedians_multiple_band_mosaic
# Select a compositing method to create your cloud-filtered mosaic
# Remove the comments from the pair of lines under one of the mosaic types
# Options are: Median or Max_NDVI
# This is the MEDIAN mosaic
scale = 0.0275
offset = 0.2
land_and_water_composite = to_f32(create_median_mosaic(land_and_water_dataset, cloud_mask), scale=scale, offset=offset)
land_composite = to_f32(create_median_mosaic(land_dataset, land_mask), scale=scale, offset=offset)
cloud_mask_composite = cloud_mask.max('time')
# This is the MAX_NDVI mosaic
# land_and_water_composite = create_max_ndvi_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_max_ndvi_mosaic(land_dataset, land_mask)
# Show the land-only composite (water will be removed ... black pixels)
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green
rgb(land_composite, bands=['swir2', 'nir', 'green'], size=8)
plt.axis('off')
plt.show()
def NDBI(dataset):
return (dataset.swir1 - dataset.nir)/(dataset.swir1 + dataset.nir)
def NDVI(dataset):
return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)
def NDWI(dataset):
return (dataset.green - dataset.nir)/(dataset.green + dataset.nir)
def SAVI(dataset):
return (dataset.nir - dataset.red)/(dataset.nir + dataset.red + 0.5)*1.5
def EVI(dataset):
return 2.5*(dataset.nir - dataset.red)/(dataset.nir + 6.0*dataset.red - 7.5*dataset.blue + 1.0)
# Water pixels masked (land only)
ndbi = NDBI(land_composite) # Normalized Difference Build Up (Urbanization) Index
ndvi = NDVI(land_composite) # Normalized Difference Vegetation Index
ndwi = NDWI(land_composite) # Normalized Difference Water Index
# Water pixels not masked
ndbi2 = NDBI(land_and_water_composite) # Normalized Difference Build Up (Urbanization) Index
ndvi2 = NDVI(land_and_water_composite) # Normalized Difference Vegetation Index
ndwi2 = NDWI(land_and_water_composite) # Normalized Difference Water Index
# Other Indices
savi = SAVI(land_composite) # Soil Adjusted Vegetation Index
evi = EVI(land_composite) # Enhanced Vegetation Index
ds_ndvi = ndvi2.to_dataset(name = "NDVI")
ds_ndwi = ndwi2.to_dataset(name= "NDWI")
ds_ndbi = ndbi2.to_dataset(name = "NDBI")
normalization_dataset = ds_ndvi.merge(ds_ndwi).merge(ds_ndbi)
Fractional Cover (FC) is used for landcover type estimation (vegetation, non-green vegetation, bare soil) of each pixel. We use a model from CSIRO (Juan Gerschmann) and apply it to a median mosaic where: Bare Soil = bs, Photosynthetic Vegetation = pv, and Non-Photosynthetic Vegetation = npv. The product is a False Color RGB result where RGB = bs/pv/npv.
land_composite = land_composite.rename_dims({'x':'longitude','y':'latitude'})
land_composite
<xarray.Dataset>
Dimensions: (latitude: 342, longitude: 322)
Coordinates:
* y (latitude) float64 4.418e+06 4.418e+06 ... 4.408e+06 4.408e+06
* x (longitude) float64 -7.386e+06 -7.386e+06 ... -7.376e+06
spatial_ref int32 6933
Dimensions without coordinates: latitude, longitude
Data variables:
red (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
green (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
blue (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
nir (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
swir1 (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
swir2 (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
pixel_qa (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>from ceos_utils.data_cube_utilities.dc_fractional_coverage_classifier import frac_coverage_classify
frac_classes = frac_coverage_classify(land_composite, clean_mask =
np.ones(land_composite.pixel_qa.shape).astype(np.bool))
/tmp/ipykernel_472/2442966445.py:3: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations np.ones(land_composite.pixel_qa.shape).astype(np.bool)) /env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. _reproject(
frac_classes
<xarray.Dataset>
Dimensions: (latitude: 342, longitude: 322)
Coordinates:
* latitude (latitude) int64 0 1 2 3 4 5 6 7 ... 335 336 337 338 339 340 341
* longitude (longitude) int64 0 1 2 3 4 5 6 7 ... 315 316 317 318 319 320 321
Data variables:
bs (latitude, longitude) float32 67.0 60.0 55.0 ... 62.0 63.0 77.0
pv (latitude, longitude) float32 22.0 28.0 32.0 ... 29.0 27.0 14.0
npv (latitude, longitude) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0# Plot of Fractional Cover
# RED = Bare Soil or Urban Areas
# BLUE = Non-Green Vegetation
# GREEN = Green Vegetation
# BLACK = Water
# Plot of RGB = NDBI-NDVI-NDWI
# RED = Bare Soil or Urban Areas
# GREEN = Vegetation
# BLUE = Water
fig, ax = plt.subplots(1, 2, figsize=(16, 8))
fc_rgb = frac_classes[['bs', 'pv', 'npv']].to_array()
ndi_rgb = normalization_dataset[['NDBI','NDVI','NDWI']].to_array()
(fc_rgb).plot.imshow(ax=ax[0], vmin=0.0, vmax=100.0)
(ndi_rgb).plot.imshow(ax=ax[1], vmin=-1.0, vmax=1.0)
ax[0].set_title('Fractional Cover'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('Normalized Difference RGB'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. _reproject(
import matplotlib as mpl
# %pylab inline
# matplotlib.rcParams['figure.figsize'] = (10, 10)
# Create a custom colour map for NDVI
# Water (blue) = NDVI -1.0 to 0.05
# Urban or Bare Soil (brown) = NDVI 0.05 to 0.25
# Low Vegetation (tan) = NDVI 0.25 to 0.4
# Croplands (light green) = NDVI 0.4 to 0.6
# Dense Vegetation / Forests (dark green) = NDVI 0.6 to 1.0
ndvi_cmap = mpl.colors.ListedColormap(['blue', '#a52a2a','#ffffcc' , '#2eb82e', '#006600'])
ndvi_bounds = [-1, 0.05, 0.25, 0.4, 0.6, 1]
ndvi_norm = mpl.colors.BoundaryNorm(ndvi_bounds, ndvi_cmap.N)
fig, ax = plt.subplots(1, 2, figsize=(16, 8))
ax[0].imshow(ndvi2, cmap="Greens", vmin=0.0, vmax=1.0)
ax[1].imshow(ndvi2, cmap=ndvi_cmap, norm = ndvi_norm)
ax[0].set_title('NDVI Basic'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('NDVI Custom Colormap'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. _reproject(
# EVI and SAVI indices using "land only" pixels
fig, ax = plt.subplots(1, 2, figsize=(18, 8))
(evi).plot.imshow(ax=ax[0], cmap="Greens", vmin=-1.0, vmax=2.5)
(savi).plot.imshow(ax=ax[1], cmap="Greens", vmin=-1.0, vmax=1.0)
ax[0].set_title('EVI'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('SAVI'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()
First we will define a minimum threshold and a maximum threshold. Then you will create a plot that colors the region between the threshold a single color (e.g. red) and the region outside the threshold will be BLACK or WHITE. Also, we will calculate the % of pixels and the number of pixels in the threshold range.
from matplotlib.ticker import FuncFormatter
def threshold_plot(da, min_threshold, max_threshold, mask = None, width = 10, *args, **kwargs):
color_in = np.array([255,0,0])
color_out = np.array([0,0,0])
color_cloud = np.array([255,255,255])
array = np.zeros((*da.values.shape, 3)).astype(np.int16)
inside = np.logical_and(da.values > min_threshold, da.values < max_threshold)
outside = np.invert(inside)
masked = np.zeros(da.values.shape).astype(bool) if mask is None else mask
array[inside] = color_in
array[outside] = color_out
array[masked] = color_cloud
def figure_ratio(ds, fixed_width = 10):
width = fixed_width
height = len(ds.latitude) * (fixed_width / len(ds.longitude))
return (width, height)
fig, ax = plt.subplots(figsize = figure_ratio(da,fixed_width = width))
lat_formatter = FuncFormatter(lambda y_val, tick_pos: "{0:.3f}".format(da.latitude.values[tick_pos] ))
lon_formatter = FuncFormatter(lambda x_val, tick_pos: "{0:.3f}".format(da.longitude.values[tick_pos]))
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.title("Threshold: {} < x < {}".format(min_threshold, max_threshold))
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.imshow(array, *args, **kwargs)
plt.axis('off')
plt.show()
# Select a threshold range for your spectral variable and generate a plot
# Remove comments from the set of 3 lines for your desired variable
# Variable choices are NDBI, NDVI, EVI, FC-Bare Soil, FC-Photosynthetic Vegetation
# NDBI (Buildup Index) = -1.0 to 1.0 (full range)
# NDBI 0.0 to 0.2 is typical for urban areas
# -----------------------
# minimum_threshold = 0.0
# maximum_threshold = 0.3
# threshold_plot(ndbi, minimum_threshold, maximum_threshold, width = 8)
# NDVI (Vegetation Index) = -1.0 to 1.0
# NDVI < 0.0 = non-vegetation (bare soil)
# NDVI 0.2 to 0.6 = grasslands
# NDVI 0.6 to 0.9 = dense vegetation / trees
# -----------------------
minimum_threshold = 0.05
maximum_threshold = 0.25
threshold_plot(ndvi.rename({'x':'longitude','y':'latitude'}), minimum_threshold, maximum_threshold, width = 8)
# EVI (Vegetation Index) = -1.0 to 2.5
# EVI 2.0 to 2.5 is typical for dense vegetation
# -----------------------
# minimum_threshold = 2.0
# maximum_threshold = 2.5
# threshold_plot(evi, minimum_threshold, maximum_threshold, width = 8)
# Fractional Cover (pv,npv,bs) = 0 to 100
# Bare Soil (bs) >40 = urbanization / bare soil
# ----------------------
# minimum_threshold = 40.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.bs, minimum_threshold, maximum_threshold, width = 8)
# Fractional Cover (pv,npv,bs) = 0 to 100
# Vegetation (pv) >80 = dense green vegetation
# ----------------------
# minimum_threshold = 80.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.pv, minimum_threshold, maximum_threshold, width = 8)
def threshold_count(da, min_threshold, max_threshold, mask = None):
def count_not_nans(arr):
return np.count_nonzero(~np.isnan(arr))
in_threshold = np.logical_and( da.values > min_threshold, da.values < max_threshold)
total_non_cloudy = count_not_nans(da.values) if mask is None else np.sum(mask.values)
return dict(total = np.size(da.values),
total_non_cloudy = total_non_cloudy,
inside = np.nansum(in_threshold),
outside = total_non_cloudy - np.nansum(in_threshold)
)
def threshold_percentage(da, min_threshold, max_threshold, mask = None):
counts = threshold_count(da, min_threshold, max_threshold, mask = mask)
return dict(percent_inside_threshold = (counts["inside"] / counts["total"]) * 100.0,
percent_outside_threshold = (counts["outside"] / counts["total"]) * 100.0,
percent_clouds = ( 100.0-counts["total_non_cloudy"] / counts["total"] * 100.0))
# Select a threshold statistical function that matches your land spectral variable
# COUNT = number of pixels in each category
# PERCENTAGE = percent of pixels in each category
# ---------------------------------
# NDBI Threshold
threshold_count(ndbi,minimum_threshold,maximum_threshold, cloud_mask_composite)
# threshold_percentage(ndbi,minimum_threshold,maximum_threshold)
# NDVI Threshold
# threshold_count(ndvi,minimum_threshold,maximum_threshold)
# threshold_percentage(ndvi,minimum_threshold,maximum_threshold)
# EVI Threshold
# threshold_count(evi,minimum_threshold,maximum_threshold)
# threshold_percentage(evi,minimum_threshold,maximum_threshold)
# Fractional Cover - Bare Soil
# threshold_count(frac_classes.bs, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.bs, minimum_threshold, maximum_threshold)
# Fractional Cover - Photosynthetic Vegetation
# threshold_count(frac_classes.pv, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.pv, minimum_threshold, maximum_threshold)
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. _reproject(
{'total': 110124, 'total_non_cloudy': 110124, 'inside': 306, 'outside': 109818}
threshold_percentage(ndbi,minimum_threshold,maximum_threshold)
{'percent_inside_threshold': 0.2778685845047401,
'percent_outside_threshold': 99.72213141549527,
'percent_clouds': 0.0}
from ceos_utils.data_cube_utilities.dc_utilities import write_geotiff_from_xr
# Remove the comment to create a GeoTIFF output product
# Change the name of the output file, or it will be overwritten for each run
# Change the desired bands at the end of the function
# Fractional Coverage
# write_geotiff_from_xr("geotiffs/frac_classes.tif", frac_classes, bands=['bs'])
# NDVI
# write_geotiff_from_xr("geotiffs/ndvi_land.tif", ndvi)
# EVI
# write_geotiff_from_xr("geotiffs/evi.tif", evi)
# WOFS
# write_geotiff_from_xr("geotiffs/wofs.tif", water_classification.wofs)
# !ls -lah geotiffs/*.tif